Bayesian GLMM Part2
Murray Logan
29 June, 2022
Preparations
Load the necessary libraries
library(rstanarm) #for fitting models in STAN
library(brms) #for fitting models in STAN
library(coda) #for diagnostics
library(bayesplot) #for diagnostics
library(ggmcmc) #for MCMC diagnostics
library(rstan) #for interfacing with STAN
library(emmeans) #for marginal means etc
library(broom) #for tidying outputs
library(DHARMa) #for residual diagnostics
library(tidybayes) #for more tidying outputs
library(ggeffects) #for partial plots
library(broom.mixed)#for tidying MCMC outputs
library(tidyverse) #for data wrangling etc
library(patchwork) #for multiple plots
library(ggridges) #for ridge plots
source('helperFunctions.R')Scenario
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
Barramundi
Format of norin.csv data files
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
Read in the data
norin <- read_csv('../public/data/norin.csv', trim_ws=TRUE)## Rows: 180 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): TRIAL
## dbl (4): FISHID, MASS, SMR_contr, CHANGE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.glimpse(norin)## Rows: 180
## Columns: 5
## $ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…We should start by declaring any categorical variables (including an varying effects) as predictors.
Exploratory data analysis
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 \sim{} \mathcal{N}(16, 35)\\ \beta_{1-6} \sim{} \mathcal{N}(0, 70)\\ \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
Lets begin by declaring the categorical predictors and random effects as factors.
norin <- norin %>% mutate(FISHID=factor(FISHID),
TRIAL=factor(TRIAL))Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.
ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
Conclusions:
- these seem reasonable enough
The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.
Doing so would introduce two additional assumptions:
- linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
- homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth(method='lm') + geom_point()## `geom_smooth()` using formula 'y ~ x'
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
geom_smooth() + geom_point()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
geom_point() + geom_line()
Conclusions:
- generally, the difference between control and trial SMR is greatest for the High temperature trial.
- there is less consistency about the relative effects of Hypoxia and Low Salinity.
- hence a random intercept/slope model might be useful.
Finally, as an acknowledgement that the metabolic response might be influenced
by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.
- add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardised by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
- add them as a regular covariate. This will include the covariate in the model and estimate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.
Lets explore the relationship between the response and Mass separately for each Trial.
ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) +
geom_point() +
geom_smooth(method='lm')## `geom_smooth()` using formula 'y ~ x'
Conclusions:
- there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
- there is no evidence that the response and Mass are related in a non-linear manner.
Fit the model
rstanarm
Using default priors
In rstanarm, the default priors are designed to be
weakly informative. They are chosen to provide moderate regularlization
(to help prevent overfitting) and help stabalise the computations.
norin.rstanarm <- stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
data = norin,
family = gaussian(),
iter = 5000,
warmup = 2000,
chains = 3,
thin = 5,
refresh = 0)norin.rstanarm %>% prior_summary()## Priors for model '.'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 20, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 20, scale = 85)
##
## Coefficients
## Specified prior:
## ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [0,0,0,...], scale = [178.99,178.99,135.15,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.03)
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more detailsThis tells us:
- for the intercept (when the family is Gaussian), a normal prior with
a mean of 20 and a standard deviation of 85 is used. A value of 2.5 is
used for scaling all parameter standard deviations. The value of 20 is
based on the mean of the response (
mean(norin$CHANGE)) and the scaled standard deviation of 85 is based on multiplying the scaling factor by the standard deviation of the response.
2.5*sd(norin$CHANGE)## [1] 84.6135- for the coefficients (in this case, just the difference between strong and weak innoculation), the default prior is a normal prior centered around 0 with a standard deviations of 178.99, 178.99, 135.15, 10.60, 34.2 and 34.2 repectively. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the stanard devation of the numerical dummy variables for the predictor (then rounded).
2.5*sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)## TRIALHypoxia TRIALLowSalinity
## 178.99305 178.99305
## SMR_contr MASS
## 135.14517 10.59372
## TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr
## 34.19156 34.19156- the auxillary prior represents whatever additional prior is required
for the nominated model. In the case of a Gaussian model, this will be
\(\sigma\), for negative binomial, it
will be the reciprocal of dispersion, for gamma, it will be shape, etc .
By default in
rstanarm, this is a exponential with a rate of 1 which is then adjusted by devision with the standard deviation of the response.
1/sd(norin$CHANGE)## [1] 0.02954611Assessing priors
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
norin.rstanarm1 <- update(norin.rstanarm, prior_PD=TRUE)norin.rstanarm1 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data=TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
#OR
norin.rstanarm1 %>%
ggemmeans(~SMR_contr*TRIAL) %>%
plot(add.data=TRUE)
Conclusions:
- we see that the range of predictions is faily wide and the predicted means could range from a small negative number to a relatively large positive number.
Defining priors
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
- \(\beta_0\): normal centered at 17
with a standard deviation of 35
- mean of 17: since
median(norin$CHANGE) - sd of 35: since
mad(norin$CHANGE)
- mean of 17: since
- \(\beta_1\): normal centred at 0
with a standard deviation of 70
- sd of 70: since
sd(norin$CHANGE) / apply(model.matrix(~scale(SMR_contr)*TRIAL + scale(MASS), norin), 2, sd)
- sd of 70: since
- \(\sigma\): exponential with rate
of 0.15
- since
1 / sd(norin$CHANGE)
- since
- \(\Sigma\): decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.
I will also overlay the raw data for comparison.
norin.rstanarm2 <- stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*scale(SMR_contr)+scale(MASS),
data = norin,
family = gaussian(),
prior_intercept = normal(17, 35, autoscale = FALSE),
prior = normal(0, 70, autoscale = FALSE),
prior_aux=rstanarm::exponential(0.03, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
prior_PD = TRUE,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)norin.rstanarm2 %>%
ggpredict(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
Now lets refit, conditioning on the data.
norin.rstanarm3 <- update(norin.rstanarm2, prior_PD=FALSE)Plotting prior and posterior
get_variables(norin.rstanarm3)## [1] "(Intercept)"
## [2] "TRIALHypoxia"
## [3] "TRIALLowSalinity"
## [4] "scale(SMR_contr)"
## [5] "scale(MASS)"
## [6] "TRIALHypoxia:scale(SMR_contr)"
## [7] "TRIALLowSalinity:scale(SMR_contr)"
## [8] "b[(Intercept) FISHID:1]"
## [9] "b[(Intercept) FISHID:2]"
## [10] "b[(Intercept) FISHID:3]"
## [11] "b[(Intercept) FISHID:4]"
## [12] "b[(Intercept) FISHID:5]"
## [13] "b[(Intercept) FISHID:6]"
## [14] "b[(Intercept) FISHID:7]"
## [15] "b[(Intercept) FISHID:8]"
## [16] "b[(Intercept) FISHID:9]"
## [17] "b[(Intercept) FISHID:10]"
## [18] "b[(Intercept) FISHID:11]"
## [19] "b[(Intercept) FISHID:12]"
## [20] "b[(Intercept) FISHID:13]"
## [21] "b[(Intercept) FISHID:14]"
## [22] "b[(Intercept) FISHID:15]"
## [23] "b[(Intercept) FISHID:16]"
## [24] "b[(Intercept) FISHID:17]"
## [25] "b[(Intercept) FISHID:18]"
## [26] "b[(Intercept) FISHID:19]"
## [27] "b[(Intercept) FISHID:20]"
## [28] "b[(Intercept) FISHID:21]"
## [29] "b[(Intercept) FISHID:22]"
## [30] "b[(Intercept) FISHID:23]"
## [31] "b[(Intercept) FISHID:24]"
## [32] "b[(Intercept) FISHID:25]"
## [33] "b[(Intercept) FISHID:26]"
## [34] "b[(Intercept) FISHID:27]"
## [35] "b[(Intercept) FISHID:28]"
## [36] "b[(Intercept) FISHID:29]"
## [37] "b[(Intercept) FISHID:30]"
## [38] "b[(Intercept) FISHID:31]"
## [39] "b[(Intercept) FISHID:32]"
## [40] "b[(Intercept) FISHID:33]"
## [41] "b[(Intercept) FISHID:34]"
## [42] "b[(Intercept) FISHID:35]"
## [43] "b[(Intercept) FISHID:36]"
## [44] "b[(Intercept) FISHID:37]"
## [45] "b[(Intercept) FISHID:38]"
## [46] "b[(Intercept) FISHID:39]"
## [47] "b[(Intercept) FISHID:40]"
## [48] "b[(Intercept) FISHID:41]"
## [49] "b[(Intercept) FISHID:42]"
## [50] "b[(Intercept) FISHID:43]"
## [51] "b[(Intercept) FISHID:44]"
## [52] "b[(Intercept) FISHID:45]"
## [53] "b[(Intercept) FISHID:46]"
## [54] "b[(Intercept) FISHID:47]"
## [55] "b[(Intercept) FISHID:48]"
## [56] "b[(Intercept) FISHID:49]"
## [57] "b[(Intercept) FISHID:50]"
## [58] "b[(Intercept) FISHID:51]"
## [59] "b[(Intercept) FISHID:52]"
## [60] "b[(Intercept) FISHID:53]"
## [61] "b[(Intercept) FISHID:54]"
## [62] "b[(Intercept) FISHID:55]"
## [63] "b[(Intercept) FISHID:56]"
## [64] "b[(Intercept) FISHID:57]"
## [65] "b[(Intercept) FISHID:58]"
## [66] "b[(Intercept) FISHID:59]"
## [67] "b[(Intercept) FISHID:60]"
## [68] "sigma"
## [69] "Sigma[FISHID:(Intercept),(Intercept)]"
## [70] "accept_stat__"
## [71] "stepsize__"
## [72] "treedepth__"
## [73] "n_leapfrog__"
## [74] "divergent__"
## [75] "energy__"posterior_vs_prior(norin.rstanarm3, color_by='vs', group_by=TRUE,
facet_args=list(scales='free_y'),
regex_pars = "^.Intercept|TRIAL|SMR_contr|MASS|sigma|Sigma")##
## Drawing from prior...
Conclusions:
- in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.
norin.rstanarm3 %>%
ggpredict(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
##OR
norin.rstanarm3 %>%
ggemmeans(~SMR_contr * TRIAL) %>%
plot(add.data = TRUE)
brms
Using default priors
In brms, the default priors are designed to be weakly
informative. They are chosen to provide moderate regularlization (to
help prevent overfitting) and help stabalise the computations.
Unlike rstanarm, brms models must be
compiled before they start sampling. For most models, the compilation of
the stan code takes around 45 seconds.
norin.form <- bf(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
family = gaussian()
)
options(width=100)
norin.form %>% get_prior(data=norin)## prior class coef group resp dpar nlpar bound
## (flat) b
## (flat) b MASS
## (flat) b SMR_contr
## (flat) b TRIALHypoxia
## (flat) b TRIALHypoxia:SMR_contr
## (flat) b TRIALLowSalinity
## (flat) b TRIALLowSalinity:SMR_contr
## student_t(3, 16.4, 34.5) Intercept
## student_t(3, 0, 34.5) sd
## student_t(3, 0, 34.5) sd FISHID
## student_t(3, 0, 34.5) sd Intercept FISHID
## student_t(3, 0, 34.5) sigma
## source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## defaultnorin.brm <- brm(norin.form,
data=norin,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0)## Compiling Stan program...## Start samplingThis tells us:
- for the intercept (when the family is Gaussian), a t-distribution
prior with a mean of 16.4 and a standard deviation of 34.5 is used:
- mean of 16.4, since
median(norin$CHANGE) - sd of 34.5 since
mad(norin$CHANGE)
- mean of 16.4, since
brmsuses flat (inproper) priors for all population effects- the prior on sigma is a half t-distribution with a mean of 0 and
standard deviation of 34.5:
- sd of 34.5 since
mad(norin$CHANGE)
- sd of 34.5 since
- the hyperprior on the standard deviation is a half t-distribution with a mean of 0 and a standard deviation of 34.5.
Defining priors
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
- \(\beta_0\): normal centred at 16
with a standard deviation of 35
- mean of 16: since
median(norin$CHANGE) - sd of 35: since
mad(norin$CHANGE)
- mean of 16: since
- \(\beta_{1-2}\) (TRIAL): normal
centred at 0 with a standard deviation of 70
- sd of 70: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 70: since
- \(\beta_{3}\) (SMR_contr): normal
centred at 0 with a standard deviation of 54
- sd of 54: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 54: since
- \(\beta_{4}\) (MASS): normal
centred at 0 with a standard deviation of 4
- sd of 4: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 4: since
- \(\beta_{5-6}\) (interaction):
normal centred at 0 with a standard deviation of 14
- sd of 14: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 14: since
- \(\sigma\): half-cauchy with
parameters 0 and 6.5
- since
sd(norin$CHANGE)
- since
- \(\sigma_j\): half-cauchy with
parameters 0 and 2.5
- since
sqrt(sd(norin$CHANGE)) - we want this prior to have most mass close to zero for the purpose of regularisation
- since
- \(\Sigma\): decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.
Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:
- half-t: as the degrees of freedom approach infinity, this will approach a half-normal
- half-cauchy: this is essentially a half-t with 1 degree of freedom
- exponential
I will also overlay the raw data for comparison.
norin %>%
group_by(TRIAL) %>%
summarise(median(CHANGE),
mad(CHANGE))norin %>%
group_by(TRIAL, FISHID) %>%
summarise(median = median(CHANGE),
MAD = mad(CHANGE)) %>%
ungroup(FISHID) %>%
summarise(sd(median))## `summarise()` has grouped output by 'TRIAL'. You can override using the
## `.groups` argument.sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)## TRIALHypoxia TRIALLowSalinity
## 71.597218 71.597218
## SMR_contr MASS
## 54.058070 4.237486
## TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr
## 13.676626 13.676626standist::visualize("normal(16,35)", xlim=c(-10,100))
standist::visualize("normal(0, 70)", "normal(0, 54)", xlim=c(-200,200))
standist::visualize("gamma(2, 1)", "gamma(35, 1)",
xlim=c(-10,100))
priors <- prior(normal(16, 35), class = 'Intercept') +
prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
prior(normal(0, 4), class = 'b', coef = 'MASS') +
prior(normal(0, 13), class = 'b') +
prior(gamma(35, 1), class = 'sigma') +
prior(cauchy(0, 2.5), class = 'sd')
norin.form <- bf(CHANGE ~ (1|FISHID) + TRIAL*SMR_contr + MASS,
family = gaussian()
)
norin.brm2 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = 'only',
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)## Compiling Stan program...## Start samplingNote in the above model, the output may have included a warning message alerting us the presence of divergent transitions. Divergent transitions are an indication that the sampler has encountered poor sampling conditions - the more divergent transitions, the more severe the issue.
Typically, divergent transitions are the result of either:
- a miss-specified model
- priors that permit the sampler to drift into unsupported areas
- complex posterior “features” (with high degrees of curvature) for which the sampler was inadequately tuned during the warmup phase
Accordingly, these divergent transitions can be addressed by either:
- reviewing the model structure
- adopting tighter priors
- increase the adaptive delta from the default of 0.8 to closer to 1. The adaptive delta defines the average acceptance probability that the sampler should aspire to during the warmup phase. Increasing the adaptive delta results in a smaller step size (and thus fewer divergences and more robust samples) however it will also result in slower sampling speeds.
norin.brm2 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
norin.brm3 <- update(norin.brm2,
sample_prior = 'yes',
control = list(adapt_delta = 0.99),
refresh = 0)## The desired updates require recompiling the model## Compiling Stan program...## Start samplingsave(norin.brm3, file = '../ws/testing/norin.brm3')norin.brm3 %>%
ggpredict(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
Plotting prior and posterior
norin.brm3 %>% get_variables()## [1] "b_Intercept" "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity" "b_SMR_contr"
## [5] "b_MASS" "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr" "sd_FISHID__Intercept"
## [9] "sigma" "r_FISHID[1,Intercept]"
## [11] "r_FISHID[2,Intercept]" "r_FISHID[3,Intercept]"
## [13] "r_FISHID[4,Intercept]" "r_FISHID[5,Intercept]"
## [15] "r_FISHID[6,Intercept]" "r_FISHID[7,Intercept]"
## [17] "r_FISHID[8,Intercept]" "r_FISHID[9,Intercept]"
## [19] "r_FISHID[10,Intercept]" "r_FISHID[11,Intercept]"
## [21] "r_FISHID[12,Intercept]" "r_FISHID[13,Intercept]"
## [23] "r_FISHID[14,Intercept]" "r_FISHID[15,Intercept]"
## [25] "r_FISHID[16,Intercept]" "r_FISHID[17,Intercept]"
## [27] "r_FISHID[18,Intercept]" "r_FISHID[19,Intercept]"
## [29] "r_FISHID[20,Intercept]" "r_FISHID[21,Intercept]"
## [31] "r_FISHID[22,Intercept]" "r_FISHID[23,Intercept]"
## [33] "r_FISHID[24,Intercept]" "r_FISHID[25,Intercept]"
## [35] "r_FISHID[26,Intercept]" "r_FISHID[27,Intercept]"
## [37] "r_FISHID[28,Intercept]" "r_FISHID[29,Intercept]"
## [39] "r_FISHID[30,Intercept]" "r_FISHID[31,Intercept]"
## [41] "r_FISHID[32,Intercept]" "r_FISHID[33,Intercept]"
## [43] "r_FISHID[34,Intercept]" "r_FISHID[35,Intercept]"
## [45] "r_FISHID[36,Intercept]" "r_FISHID[37,Intercept]"
## [47] "r_FISHID[38,Intercept]" "r_FISHID[39,Intercept]"
## [49] "r_FISHID[40,Intercept]" "r_FISHID[41,Intercept]"
## [51] "r_FISHID[42,Intercept]" "r_FISHID[43,Intercept]"
## [53] "r_FISHID[44,Intercept]" "r_FISHID[45,Intercept]"
## [55] "r_FISHID[46,Intercept]" "r_FISHID[47,Intercept]"
## [57] "r_FISHID[48,Intercept]" "r_FISHID[49,Intercept]"
## [59] "r_FISHID[50,Intercept]" "r_FISHID[51,Intercept]"
## [61] "r_FISHID[52,Intercept]" "r_FISHID[53,Intercept]"
## [63] "r_FISHID[54,Intercept]" "r_FISHID[55,Intercept]"
## [65] "r_FISHID[56,Intercept]" "r_FISHID[57,Intercept]"
## [67] "r_FISHID[58,Intercept]" "r_FISHID[59,Intercept]"
## [69] "r_FISHID[60,Intercept]" "prior_Intercept"
## [71] "prior_b_TRIALHypoxia" "prior_b_TRIALLowSalinity"
## [73] "prior_b_SMR_contr" "prior_b_MASS"
## [75] "prior_b_TRIALHypoxia:SMR_contr" "prior_b_TRIALLowSalinity:SMR_contr"
## [77] "prior_sigma" "prior_sd_FISHID"
## [79] "lp__" "accept_stat__"
## [81] "stepsize__" "treedepth__"
## [83] "n_leapfrog__" "divergent__"
## [85] "energy__"norin.brm3 %>% hypothesis('TRIALHypoxia=0') %>% plot
norin.brm3 %>% hypothesis('SMR_contr=0') %>% plot
norin.brm3 %>% SUYR_prior_and_posterior()## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
Random intercept/slope model
While we are here, we might like to explore a random intercept/slope model
priors <- prior(normal(16, 35), class = 'Intercept') +
prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
prior(normal(0, 4), class = 'b', coef = 'MASS') +
prior(normal(0, 13), class = 'b') +
prior(gamma(35, 1), class = 'sigma') +
prior(cauchy(0, 2.5), class = 'sd') +
prior(lkj_corr_cholesky(1), class = 'L')
norin.form <- bf(CHANGE ~ (TRIAL|FISHID) + TRIAL*SMR_contr + MASS,
family = gaussian()
)
norin.brm4 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = 'yes',
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0,
control = list(adapt_delta=0.99)
)## Compiling Stan program...## Start samplingsave(norin.brm4, file = '../ws/testing/norin.brm4')norin.brm4 %>% get_variables()## [1] "b_Intercept"
## [2] "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity"
## [4] "b_SMR_contr"
## [5] "b_MASS"
## [6] "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr"
## [8] "sd_FISHID__Intercept"
## [9] "sd_FISHID__TRIALHypoxia"
## [10] "sd_FISHID__TRIALLowSalinity"
## [11] "cor_FISHID__Intercept__TRIALHypoxia"
## [12] "cor_FISHID__Intercept__TRIALLowSalinity"
## [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
## [14] "sigma"
## [15] "r_FISHID[1,Intercept]"
## [16] "r_FISHID[2,Intercept]"
## [17] "r_FISHID[3,Intercept]"
## [18] "r_FISHID[4,Intercept]"
## [19] "r_FISHID[5,Intercept]"
## [20] "r_FISHID[6,Intercept]"
## [21] "r_FISHID[7,Intercept]"
## [22] "r_FISHID[8,Intercept]"
## [23] "r_FISHID[9,Intercept]"
## [24] "r_FISHID[10,Intercept]"
## [25] "r_FISHID[11,Intercept]"
## [26] "r_FISHID[12,Intercept]"
## [27] "r_FISHID[13,Intercept]"
## [28] "r_FISHID[14,Intercept]"
## [29] "r_FISHID[15,Intercept]"
## [30] "r_FISHID[16,Intercept]"
## [31] "r_FISHID[17,Intercept]"
## [32] "r_FISHID[18,Intercept]"
## [33] "r_FISHID[19,Intercept]"
## [34] "r_FISHID[20,Intercept]"
## [35] "r_FISHID[21,Intercept]"
## [36] "r_FISHID[22,Intercept]"
## [37] "r_FISHID[23,Intercept]"
## [38] "r_FISHID[24,Intercept]"
## [39] "r_FISHID[25,Intercept]"
## [40] "r_FISHID[26,Intercept]"
## [41] "r_FISHID[27,Intercept]"
## [42] "r_FISHID[28,Intercept]"
## [43] "r_FISHID[29,Intercept]"
## [44] "r_FISHID[30,Intercept]"
## [45] "r_FISHID[31,Intercept]"
## [46] "r_FISHID[32,Intercept]"
## [47] "r_FISHID[33,Intercept]"
## [48] "r_FISHID[34,Intercept]"
## [49] "r_FISHID[35,Intercept]"
## [50] "r_FISHID[36,Intercept]"
## [51] "r_FISHID[37,Intercept]"
## [52] "r_FISHID[38,Intercept]"
## [53] "r_FISHID[39,Intercept]"
## [54] "r_FISHID[40,Intercept]"
## [55] "r_FISHID[41,Intercept]"
## [56] "r_FISHID[42,Intercept]"
## [57] "r_FISHID[43,Intercept]"
## [58] "r_FISHID[44,Intercept]"
## [59] "r_FISHID[45,Intercept]"
## [60] "r_FISHID[46,Intercept]"
## [61] "r_FISHID[47,Intercept]"
## [62] "r_FISHID[48,Intercept]"
## [63] "r_FISHID[49,Intercept]"
## [64] "r_FISHID[50,Intercept]"
## [65] "r_FISHID[51,Intercept]"
## [66] "r_FISHID[52,Intercept]"
## [67] "r_FISHID[53,Intercept]"
## [68] "r_FISHID[54,Intercept]"
## [69] "r_FISHID[55,Intercept]"
## [70] "r_FISHID[56,Intercept]"
## [71] "r_FISHID[57,Intercept]"
## [72] "r_FISHID[58,Intercept]"
## [73] "r_FISHID[59,Intercept]"
## [74] "r_FISHID[60,Intercept]"
## [75] "r_FISHID[1,TRIALHypoxia]"
## [76] "r_FISHID[2,TRIALHypoxia]"
## [77] "r_FISHID[3,TRIALHypoxia]"
## [78] "r_FISHID[4,TRIALHypoxia]"
## [79] "r_FISHID[5,TRIALHypoxia]"
## [80] "r_FISHID[6,TRIALHypoxia]"
## [81] "r_FISHID[7,TRIALHypoxia]"
## [82] "r_FISHID[8,TRIALHypoxia]"
## [83] "r_FISHID[9,TRIALHypoxia]"
## [84] "r_FISHID[10,TRIALHypoxia]"
## [85] "r_FISHID[11,TRIALHypoxia]"
## [86] "r_FISHID[12,TRIALHypoxia]"
## [87] "r_FISHID[13,TRIALHypoxia]"
## [88] "r_FISHID[14,TRIALHypoxia]"
## [89] "r_FISHID[15,TRIALHypoxia]"
## [90] "r_FISHID[16,TRIALHypoxia]"
## [91] "r_FISHID[17,TRIALHypoxia]"
## [92] "r_FISHID[18,TRIALHypoxia]"
## [93] "r_FISHID[19,TRIALHypoxia]"
## [94] "r_FISHID[20,TRIALHypoxia]"
## [95] "r_FISHID[21,TRIALHypoxia]"
## [96] "r_FISHID[22,TRIALHypoxia]"
## [97] "r_FISHID[23,TRIALHypoxia]"
## [98] "r_FISHID[24,TRIALHypoxia]"
## [99] "r_FISHID[25,TRIALHypoxia]"
## [100] "r_FISHID[26,TRIALHypoxia]"
## [101] "r_FISHID[27,TRIALHypoxia]"
## [102] "r_FISHID[28,TRIALHypoxia]"
## [103] "r_FISHID[29,TRIALHypoxia]"
## [104] "r_FISHID[30,TRIALHypoxia]"
## [105] "r_FISHID[31,TRIALHypoxia]"
## [106] "r_FISHID[32,TRIALHypoxia]"
## [107] "r_FISHID[33,TRIALHypoxia]"
## [108] "r_FISHID[34,TRIALHypoxia]"
## [109] "r_FISHID[35,TRIALHypoxia]"
## [110] "r_FISHID[36,TRIALHypoxia]"
## [111] "r_FISHID[37,TRIALHypoxia]"
## [112] "r_FISHID[38,TRIALHypoxia]"
## [113] "r_FISHID[39,TRIALHypoxia]"
## [114] "r_FISHID[40,TRIALHypoxia]"
## [115] "r_FISHID[41,TRIALHypoxia]"
## [116] "r_FISHID[42,TRIALHypoxia]"
## [117] "r_FISHID[43,TRIALHypoxia]"
## [118] "r_FISHID[44,TRIALHypoxia]"
## [119] "r_FISHID[45,TRIALHypoxia]"
## [120] "r_FISHID[46,TRIALHypoxia]"
## [121] "r_FISHID[47,TRIALHypoxia]"
## [122] "r_FISHID[48,TRIALHypoxia]"
## [123] "r_FISHID[49,TRIALHypoxia]"
## [124] "r_FISHID[50,TRIALHypoxia]"
## [125] "r_FISHID[51,TRIALHypoxia]"
## [126] "r_FISHID[52,TRIALHypoxia]"
## [127] "r_FISHID[53,TRIALHypoxia]"
## [128] "r_FISHID[54,TRIALHypoxia]"
## [129] "r_FISHID[55,TRIALHypoxia]"
## [130] "r_FISHID[56,TRIALHypoxia]"
## [131] "r_FISHID[57,TRIALHypoxia]"
## [132] "r_FISHID[58,TRIALHypoxia]"
## [133] "r_FISHID[59,TRIALHypoxia]"
## [134] "r_FISHID[60,TRIALHypoxia]"
## [135] "r_FISHID[1,TRIALLowSalinity]"
## [136] "r_FISHID[2,TRIALLowSalinity]"
## [137] "r_FISHID[3,TRIALLowSalinity]"
## [138] "r_FISHID[4,TRIALLowSalinity]"
## [139] "r_FISHID[5,TRIALLowSalinity]"
## [140] "r_FISHID[6,TRIALLowSalinity]"
## [141] "r_FISHID[7,TRIALLowSalinity]"
## [142] "r_FISHID[8,TRIALLowSalinity]"
## [143] "r_FISHID[9,TRIALLowSalinity]"
## [144] "r_FISHID[10,TRIALLowSalinity]"
## [145] "r_FISHID[11,TRIALLowSalinity]"
## [146] "r_FISHID[12,TRIALLowSalinity]"
## [147] "r_FISHID[13,TRIALLowSalinity]"
## [148] "r_FISHID[14,TRIALLowSalinity]"
## [149] "r_FISHID[15,TRIALLowSalinity]"
## [150] "r_FISHID[16,TRIALLowSalinity]"
## [151] "r_FISHID[17,TRIALLowSalinity]"
## [152] "r_FISHID[18,TRIALLowSalinity]"
## [153] "r_FISHID[19,TRIALLowSalinity]"
## [154] "r_FISHID[20,TRIALLowSalinity]"
## [155] "r_FISHID[21,TRIALLowSalinity]"
## [156] "r_FISHID[22,TRIALLowSalinity]"
## [157] "r_FISHID[23,TRIALLowSalinity]"
## [158] "r_FISHID[24,TRIALLowSalinity]"
## [159] "r_FISHID[25,TRIALLowSalinity]"
## [160] "r_FISHID[26,TRIALLowSalinity]"
## [161] "r_FISHID[27,TRIALLowSalinity]"
## [162] "r_FISHID[28,TRIALLowSalinity]"
## [163] "r_FISHID[29,TRIALLowSalinity]"
## [164] "r_FISHID[30,TRIALLowSalinity]"
## [165] "r_FISHID[31,TRIALLowSalinity]"
## [166] "r_FISHID[32,TRIALLowSalinity]"
## [167] "r_FISHID[33,TRIALLowSalinity]"
## [168] "r_FISHID[34,TRIALLowSalinity]"
## [169] "r_FISHID[35,TRIALLowSalinity]"
## [170] "r_FISHID[36,TRIALLowSalinity]"
## [171] "r_FISHID[37,TRIALLowSalinity]"
## [172] "r_FISHID[38,TRIALLowSalinity]"
## [173] "r_FISHID[39,TRIALLowSalinity]"
## [174] "r_FISHID[40,TRIALLowSalinity]"
## [175] "r_FISHID[41,TRIALLowSalinity]"
## [176] "r_FISHID[42,TRIALLowSalinity]"
## [177] "r_FISHID[43,TRIALLowSalinity]"
## [178] "r_FISHID[44,TRIALLowSalinity]"
## [179] "r_FISHID[45,TRIALLowSalinity]"
## [180] "r_FISHID[46,TRIALLowSalinity]"
## [181] "r_FISHID[47,TRIALLowSalinity]"
## [182] "r_FISHID[48,TRIALLowSalinity]"
## [183] "r_FISHID[49,TRIALLowSalinity]"
## [184] "r_FISHID[50,TRIALLowSalinity]"
## [185] "r_FISHID[51,TRIALLowSalinity]"
## [186] "r_FISHID[52,TRIALLowSalinity]"
## [187] "r_FISHID[53,TRIALLowSalinity]"
## [188] "r_FISHID[54,TRIALLowSalinity]"
## [189] "r_FISHID[55,TRIALLowSalinity]"
## [190] "r_FISHID[56,TRIALLowSalinity]"
## [191] "r_FISHID[57,TRIALLowSalinity]"
## [192] "r_FISHID[58,TRIALLowSalinity]"
## [193] "r_FISHID[59,TRIALLowSalinity]"
## [194] "r_FISHID[60,TRIALLowSalinity]"
## [195] "prior_Intercept"
## [196] "prior_b_TRIALHypoxia"
## [197] "prior_b_TRIALLowSalinity"
## [198] "prior_b_SMR_contr"
## [199] "prior_b_MASS"
## [200] "prior_b_TRIALHypoxia:SMR_contr"
## [201] "prior_b_TRIALLowSalinity:SMR_contr"
## [202] "prior_sigma"
## [203] "prior_sd_FISHID"
## [204] "prior_cor_FISHID"
## [205] "lp__"
## [206] "accept_stat__"
## [207] "stepsize__"
## [208] "treedepth__"
## [209] "n_leapfrog__"
## [210] "divergent__"
## [211] "energy__"norin.brm4 %>% hypothesis('TRIALHypoxia=0') %>% plot
norin.brm4 %>% hypothesis('SMR_contr=0') %>% plot
norin.brm4 %>% SUYR_prior_and_posterior()## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
norin.brm4 %>%
posterior_samples %>%
dplyr::select(-`lp__`) %>%
pivot_longer(everything(), names_to = 'key') %>%
filter(!str_detect(key, '^r')) %>%
mutate(Type = ifelse(str_detect(key, 'prior'), 'Prior', 'Posterior'),
Class = case_when(
str_detect(key, '(^b|^prior).*Intercept$') ~ 'Intercept',
str_detect(key, 'b_TRIAL.*|prior_b_TRIAL.*') & str_detect(key, '.*\\:.*') ~ 'TRIAL',
str_detect(key, 'b_SMR_contr|prior_b_SMR_contr') ~ 'SMR',
str_detect(key, 'b_MASS|prior_b_MASS') ~ 'MASS',
str_detect(key, '.*\\:.*|prior_b_.*\\:.*') ~ 'Interaction',
str_detect(key, 'sd') ~ 'sd',
str_detect(key, '^cor|prior_cor') ~ 'cor',
str_detect(key, 'sigma') ~ 'sigma'),
Par = str_replace(key, 'b_', '')) %>%
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge())+
facet_wrap(~Class, scales = 'free')## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
Compare models
We can compare the two models using LOO
(l.1 <- norin.brm3 %>% loo())##
## Computed from 2400 by 180 log-likelihood matrix
##
## Estimate SE
## elpd_loo -833.4 12.8
## p_loo 15.2 2.7
## looic 1666.7 25.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 177 98.3% 275
## (0.5, 0.7] (ok) 3 1.7% 467
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.(l.2 <- norin.brm4 %>% loo())## Warning: Found 35 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.##
## Computed from 2400 by 180 log-likelihood matrix
##
## Estimate SE
## elpd_loo -807.0 9.5
## p_loo 54.5 6.1
## looic 1613.9 19.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 96 53.3% 522
## (0.5, 0.7] (ok) 49 27.2% 128
## (0.7, 1] (bad) 30 16.7% 21
## (1, Inf) (very bad) 5 2.8% 7
## See help('pareto-k-diagnostic') for details.loo_compare(l.1, l.2)## elpd_diff se_diff
## . 0.0 0.0
## . -26.4 6.6There is substantially more support for the more complex random intercept/slope model over the simpler random intercept only model. Consequently, we will continue with the random intercept/slope model.
MCMC sampling diagnostics
In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.
There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.
brms
bayesplot
The bayesplot package offers a range of MCMC diagnostics
as well as Posterior Probability Checks (PPC), all of which have a
convenient plot() interface. Lets start with the MCMC
diagnostics.
See list of available diagnostics by name
available_mcmc()## bayesplot MCMC module:
## mcmc_acf
## mcmc_acf_bar
## mcmc_areas
## mcmc_areas_data
## mcmc_areas_ridges
## mcmc_areas_ridges_data
## mcmc_combo
## mcmc_dens
## mcmc_dens_chains
## mcmc_dens_chains_data
## mcmc_dens_overlay
## mcmc_hex
## mcmc_hist
## mcmc_hist_by_chain
## mcmc_intervals
## mcmc_intervals_data
## mcmc_neff
## mcmc_neff_data
## mcmc_neff_hist
## mcmc_nuts_acceptance
## mcmc_nuts_divergence
## mcmc_nuts_energy
## mcmc_nuts_stepsize
## mcmc_nuts_treedepth
## mcmc_pairs
## mcmc_parcoord
## mcmc_parcoord_data
## mcmc_rank_hist
## mcmc_rank_overlay
## mcmc_recover_hist
## mcmc_recover_intervals
## mcmc_recover_scatter
## mcmc_rhat
## mcmc_rhat_data
## mcmc_rhat_hist
## mcmc_scatter
## mcmc_trace
## mcmc_trace_data
## mcmc_trace_highlight
## mcmc_violinOf these, we will focus on:
- trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
norin.brm4 %>% mcmc_plot(type='trace')## No divergences to plot.
The chains appear well mixed and very similar
- acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
norin.brm4 %>% mcmc_plot(type='acf_bar')
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
- rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
norin.brm4 %>% mcmc_plot(type='rhat_hist')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
norin.brm4 %>% mcmc_plot(type='neff_hist')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
More diagnostics
norin.brm4 %>% mcmc_plot(type='combo')
norin.brm4 %>% mcmc_plot(type='violin')
stan plots
The rstan package offers a range of MCMC diagnostics.
Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
norin.brm4 %>% get_variables()## [1] "b_Intercept"
## [2] "b_TRIALHypoxia"
## [3] "b_TRIALLowSalinity"
## [4] "b_SMR_contr"
## [5] "b_MASS"
## [6] "b_TRIALHypoxia:SMR_contr"
## [7] "b_TRIALLowSalinity:SMR_contr"
## [8] "sd_FISHID__Intercept"
## [9] "sd_FISHID__TRIALHypoxia"
## [10] "sd_FISHID__TRIALLowSalinity"
## [11] "cor_FISHID__Intercept__TRIALHypoxia"
## [12] "cor_FISHID__Intercept__TRIALLowSalinity"
## [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
## [14] "sigma"
## [15] "r_FISHID[1,Intercept]"
## [16] "r_FISHID[2,Intercept]"
## [17] "r_FISHID[3,Intercept]"
## [18] "r_FISHID[4,Intercept]"
## [19] "r_FISHID[5,Intercept]"
## [20] "r_FISHID[6,Intercept]"
## [21] "r_FISHID[7,Intercept]"
## [22] "r_FISHID[8,Intercept]"
## [23] "r_FISHID[9,Intercept]"
## [24] "r_FISHID[10,Intercept]"
## [25] "r_FISHID[11,Intercept]"
## [26] "r_FISHID[12,Intercept]"
## [27] "r_FISHID[13,Intercept]"
## [28] "r_FISHID[14,Intercept]"
## [29] "r_FISHID[15,Intercept]"
## [30] "r_FISHID[16,Intercept]"
## [31] "r_FISHID[17,Intercept]"
## [32] "r_FISHID[18,Intercept]"
## [33] "r_FISHID[19,Intercept]"
## [34] "r_FISHID[20,Intercept]"
## [35] "r_FISHID[21,Intercept]"
## [36] "r_FISHID[22,Intercept]"
## [37] "r_FISHID[23,Intercept]"
## [38] "r_FISHID[24,Intercept]"
## [39] "r_FISHID[25,Intercept]"
## [40] "r_FISHID[26,Intercept]"
## [41] "r_FISHID[27,Intercept]"
## [42] "r_FISHID[28,Intercept]"
## [43] "r_FISHID[29,Intercept]"
## [44] "r_FISHID[30,Intercept]"
## [45] "r_FISHID[31,Intercept]"
## [46] "r_FISHID[32,Intercept]"
## [47] "r_FISHID[33,Intercept]"
## [48] "r_FISHID[34,Intercept]"
## [49] "r_FISHID[35,Intercept]"
## [50] "r_FISHID[36,Intercept]"
## [51] "r_FISHID[37,Intercept]"
## [52] "r_FISHID[38,Intercept]"
## [53] "r_FISHID[39,Intercept]"
## [54] "r_FISHID[40,Intercept]"
## [55] "r_FISHID[41,Intercept]"
## [56] "r_FISHID[42,Intercept]"
## [57] "r_FISHID[43,Intercept]"
## [58] "r_FISHID[44,Intercept]"
## [59] "r_FISHID[45,Intercept]"
## [60] "r_FISHID[46,Intercept]"
## [61] "r_FISHID[47,Intercept]"
## [62] "r_FISHID[48,Intercept]"
## [63] "r_FISHID[49,Intercept]"
## [64] "r_FISHID[50,Intercept]"
## [65] "r_FISHID[51,Intercept]"
## [66] "r_FISHID[52,Intercept]"
## [67] "r_FISHID[53,Intercept]"
## [68] "r_FISHID[54,Intercept]"
## [69] "r_FISHID[55,Intercept]"
## [70] "r_FISHID[56,Intercept]"
## [71] "r_FISHID[57,Intercept]"
## [72] "r_FISHID[58,Intercept]"
## [73] "r_FISHID[59,Intercept]"
## [74] "r_FISHID[60,Intercept]"
## [75] "r_FISHID[1,TRIALHypoxia]"
## [76] "r_FISHID[2,TRIALHypoxia]"
## [77] "r_FISHID[3,TRIALHypoxia]"
## [78] "r_FISHID[4,TRIALHypoxia]"
## [79] "r_FISHID[5,TRIALHypoxia]"
## [80] "r_FISHID[6,TRIALHypoxia]"
## [81] "r_FISHID[7,TRIALHypoxia]"
## [82] "r_FISHID[8,TRIALHypoxia]"
## [83] "r_FISHID[9,TRIALHypoxia]"
## [84] "r_FISHID[10,TRIALHypoxia]"
## [85] "r_FISHID[11,TRIALHypoxia]"
## [86] "r_FISHID[12,TRIALHypoxia]"
## [87] "r_FISHID[13,TRIALHypoxia]"
## [88] "r_FISHID[14,TRIALHypoxia]"
## [89] "r_FISHID[15,TRIALHypoxia]"
## [90] "r_FISHID[16,TRIALHypoxia]"
## [91] "r_FISHID[17,TRIALHypoxia]"
## [92] "r_FISHID[18,TRIALHypoxia]"
## [93] "r_FISHID[19,TRIALHypoxia]"
## [94] "r_FISHID[20,TRIALHypoxia]"
## [95] "r_FISHID[21,TRIALHypoxia]"
## [96] "r_FISHID[22,TRIALHypoxia]"
## [97] "r_FISHID[23,TRIALHypoxia]"
## [98] "r_FISHID[24,TRIALHypoxia]"
## [99] "r_FISHID[25,TRIALHypoxia]"
## [100] "r_FISHID[26,TRIALHypoxia]"
## [101] "r_FISHID[27,TRIALHypoxia]"
## [102] "r_FISHID[28,TRIALHypoxia]"
## [103] "r_FISHID[29,TRIALHypoxia]"
## [104] "r_FISHID[30,TRIALHypoxia]"
## [105] "r_FISHID[31,TRIALHypoxia]"
## [106] "r_FISHID[32,TRIALHypoxia]"
## [107] "r_FISHID[33,TRIALHypoxia]"
## [108] "r_FISHID[34,TRIALHypoxia]"
## [109] "r_FISHID[35,TRIALHypoxia]"
## [110] "r_FISHID[36,TRIALHypoxia]"
## [111] "r_FISHID[37,TRIALHypoxia]"
## [112] "r_FISHID[38,TRIALHypoxia]"
## [113] "r_FISHID[39,TRIALHypoxia]"
## [114] "r_FISHID[40,TRIALHypoxia]"
## [115] "r_FISHID[41,TRIALHypoxia]"
## [116] "r_FISHID[42,TRIALHypoxia]"
## [117] "r_FISHID[43,TRIALHypoxia]"
## [118] "r_FISHID[44,TRIALHypoxia]"
## [119] "r_FISHID[45,TRIALHypoxia]"
## [120] "r_FISHID[46,TRIALHypoxia]"
## [121] "r_FISHID[47,TRIALHypoxia]"
## [122] "r_FISHID[48,TRIALHypoxia]"
## [123] "r_FISHID[49,TRIALHypoxia]"
## [124] "r_FISHID[50,TRIALHypoxia]"
## [125] "r_FISHID[51,TRIALHypoxia]"
## [126] "r_FISHID[52,TRIALHypoxia]"
## [127] "r_FISHID[53,TRIALHypoxia]"
## [128] "r_FISHID[54,TRIALHypoxia]"
## [129] "r_FISHID[55,TRIALHypoxia]"
## [130] "r_FISHID[56,TRIALHypoxia]"
## [131] "r_FISHID[57,TRIALHypoxia]"
## [132] "r_FISHID[58,TRIALHypoxia]"
## [133] "r_FISHID[59,TRIALHypoxia]"
## [134] "r_FISHID[60,TRIALHypoxia]"
## [135] "r_FISHID[1,TRIALLowSalinity]"
## [136] "r_FISHID[2,TRIALLowSalinity]"
## [137] "r_FISHID[3,TRIALLowSalinity]"
## [138] "r_FISHID[4,TRIALLowSalinity]"
## [139] "r_FISHID[5,TRIALLowSalinity]"
## [140] "r_FISHID[6,TRIALLowSalinity]"
## [141] "r_FISHID[7,TRIALLowSalinity]"
## [142] "r_FISHID[8,TRIALLowSalinity]"
## [143] "r_FISHID[9,TRIALLowSalinity]"
## [144] "r_FISHID[10,TRIALLowSalinity]"
## [145] "r_FISHID[11,TRIALLowSalinity]"
## [146] "r_FISHID[12,TRIALLowSalinity]"
## [147] "r_FISHID[13,TRIALLowSalinity]"
## [148] "r_FISHID[14,TRIALLowSalinity]"
## [149] "r_FISHID[15,TRIALLowSalinity]"
## [150] "r_FISHID[16,TRIALLowSalinity]"
## [151] "r_FISHID[17,TRIALLowSalinity]"
## [152] "r_FISHID[18,TRIALLowSalinity]"
## [153] "r_FISHID[19,TRIALLowSalinity]"
## [154] "r_FISHID[20,TRIALLowSalinity]"
## [155] "r_FISHID[21,TRIALLowSalinity]"
## [156] "r_FISHID[22,TRIALLowSalinity]"
## [157] "r_FISHID[23,TRIALLowSalinity]"
## [158] "r_FISHID[24,TRIALLowSalinity]"
## [159] "r_FISHID[25,TRIALLowSalinity]"
## [160] "r_FISHID[26,TRIALLowSalinity]"
## [161] "r_FISHID[27,TRIALLowSalinity]"
## [162] "r_FISHID[28,TRIALLowSalinity]"
## [163] "r_FISHID[29,TRIALLowSalinity]"
## [164] "r_FISHID[30,TRIALLowSalinity]"
## [165] "r_FISHID[31,TRIALLowSalinity]"
## [166] "r_FISHID[32,TRIALLowSalinity]"
## [167] "r_FISHID[33,TRIALLowSalinity]"
## [168] "r_FISHID[34,TRIALLowSalinity]"
## [169] "r_FISHID[35,TRIALLowSalinity]"
## [170] "r_FISHID[36,TRIALLowSalinity]"
## [171] "r_FISHID[37,TRIALLowSalinity]"
## [172] "r_FISHID[38,TRIALLowSalinity]"
## [173] "r_FISHID[39,TRIALLowSalinity]"
## [174] "r_FISHID[40,TRIALLowSalinity]"
## [175] "r_FISHID[41,TRIALLowSalinity]"
## [176] "r_FISHID[42,TRIALLowSalinity]"
## [177] "r_FISHID[43,TRIALLowSalinity]"
## [178] "r_FISHID[44,TRIALLowSalinity]"
## [179] "r_FISHID[45,TRIALLowSalinity]"
## [180] "r_FISHID[46,TRIALLowSalinity]"
## [181] "r_FISHID[47,TRIALLowSalinity]"
## [182] "r_FISHID[48,TRIALLowSalinity]"
## [183] "r_FISHID[49,TRIALLowSalinity]"
## [184] "r_FISHID[50,TRIALLowSalinity]"
## [185] "r_FISHID[51,TRIALLowSalinity]"
## [186] "r_FISHID[52,TRIALLowSalinity]"
## [187] "r_FISHID[53,TRIALLowSalinity]"
## [188] "r_FISHID[54,TRIALLowSalinity]"
## [189] "r_FISHID[55,TRIALLowSalinity]"
## [190] "r_FISHID[56,TRIALLowSalinity]"
## [191] "r_FISHID[57,TRIALLowSalinity]"
## [192] "r_FISHID[58,TRIALLowSalinity]"
## [193] "r_FISHID[59,TRIALLowSalinity]"
## [194] "r_FISHID[60,TRIALLowSalinity]"
## [195] "prior_Intercept"
## [196] "prior_b_TRIALHypoxia"
## [197] "prior_b_TRIALLowSalinity"
## [198] "prior_b_SMR_contr"
## [199] "prior_b_MASS"
## [200] "prior_b_TRIALHypoxia:SMR_contr"
## [201] "prior_b_TRIALLowSalinity:SMR_contr"
## [202] "prior_sigma"
## [203] "prior_sd_FISHID"
## [204] "prior_cor_FISHID"
## [205] "lp__"
## [206] "accept_stat__"
## [207] "stepsize__"
## [208] "treedepth__"
## [209] "n_leapfrog__"
## [210] "divergent__"
## [211] "energy__"pars <- norin.brm4 %>% get_variables()
pars <- str_extract(pars, '^b_.*|^sigma$|^sd.*') %>% na.omit()
norin.brm4$fit %>%
stan_trace(pars = pars)
The chains appear well mixed and very similar
- stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
norin.brm4$fit %>%
stan_ac(pars = pars)
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
norin.brm4$fit %>% stan_rhat() ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
norin.brm4$fit %>% stan_ess()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
norin.brm4$fit %>%
stan_dens(separate_chains = TRUE, pars = pars)
ggmcmc
The ggmean package also as a set of MCMC diagnostic
functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
norin.ggs <- norin.brm4 %>% ggs(burnin = FALSE, inc_warmup = FALSE)## Warning in custom.sort(D$Parameter): NAs introduced by coercionnorin.ggs %>% ggs_traceplot()
The chains appear well mixed and very similar
- gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
ggs_autocorrelation(norin.ggs)
There is no evidence of autocorrelation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
ggs_Rhat(norin.ggs, scaling = 1.01)
All Rhat values are below 1.01, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
ggs_effective(norin.ggs)
Ratios all very high.
More diagnostics
ggs_crosscorrelation(norin.ggs)
ggs_grb(norin.ggs)
pp check
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
See list of available diagnostics by name
available_ppc()## bayesplot PPC module:
## ppc_bars
## ppc_bars_grouped
## ppc_boxplot
## ppc_data
## ppc_dens
## ppc_dens_overlay
## ppc_dens_overlay_grouped
## ppc_ecdf_overlay
## ppc_ecdf_overlay_grouped
## ppc_error_binned
## ppc_error_hist
## ppc_error_hist_grouped
## ppc_error_scatter
## ppc_error_scatter_avg
## ppc_error_scatter_avg_vs_x
## ppc_freqpoly
## ppc_freqpoly_grouped
## ppc_hist
## ppc_intervals
## ppc_intervals_data
## ppc_intervals_grouped
## ppc_km_overlay
## ppc_loo_intervals
## ppc_loo_pit
## ppc_loo_pit_data
## ppc_loo_pit_overlay
## ppc_loo_pit_qq
## ppc_loo_ribbon
## ppc_ribbon
## ppc_ribbon_data
## ppc_ribbon_grouped
## ppc_rootogram
## ppc_scatter
## ppc_scatter_avg
## ppc_scatter_avg_grouped
## ppc_stat
## ppc_stat_2d
## ppc_stat_freqpoly_grouped
## ppc_stat_grouped
## ppc_violin_grouped- dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
norin.brm4 %>% pp_check(type = 'dens_overlay', ndraws = 100)
The model draws appear to be consistent with the observed data.
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response
norin.brm4 %>% pp_check(type = 'error_scatter_avg')## Using all posterior draws for ppc type 'error_scatter_avg' by default.
This is not really interpretable
- intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
norin.brm4 %>% pp_check(group = 'TRIAL', type = 'intervals')## Using all posterior draws for ppc type 'intervals' by default.## Warning: The following arguments were unrecognized and ignored: group
norin.brm3 %>% pp_check(group = 'TRIAL', type = 'intervals_grouped')## Using all posterior draws for ppc type 'intervals_grouped' by default.
norin.brm3 %>% pp_check(group = 'TRIAL', type = 'violin_grouped')## Using all posterior draws for ppc type 'violin_grouped' by default.
The shinystan package allows the full suite of MCMC
diagnostics and posterior predictive checks to be accessed via a web
interface.
#library(shinystan)
#launch_shinystan(norin.brm2)DHARMa residuals
DHARMa residuals provide very useful diagnostics. Unfortunately, we
cannot directly use the simulateResiduals() function to
generate the simulated residuals. However, if we are willing to
calculate some of the components ourself, we can still obtain the
simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
preds <- norin.brm4 %>% posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE)
plot(norin.resids, quantreg = FALSE)
Conclusions:
- the simulated residuals do not suggest any issues with the residuals
- there is no evidence of a lack of fit.
Partial effects plots
brms
conditional_effects
norin.brm4 %>%
conditional_effects() %>%
plot(ask = FALSE, points = TRUE, plot = FALSE) %>%
wrap_plots()
ggpredict
norin.brm4 %>%
ggpredict(terms = c("SMR_contr", "TRIAL")) %>%
plot(add.data = TRUE)## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.
ggemmeans
norin.brm4 %>%
ggemmeans(~SMR_contr*TRIAL) %>%
plot(add.data = TRUE)
fitted_draws
It is not really possible to do this via the fitted draws as it would
not be marginalising over MASS or FISHID.
Model investigation
brms
brms captures the MCMC samples from stan
within the returned list. There are numerous ways to retrieve and
summarise these samples. The first three provide convenient numeric
summaries from which you can draw conclusions, the last four provide
ways of obtaining the full posteriors.
summary
The summary() method generates simple summaries (mean,
standard deviation as well as 10, 50 and 90 percentiles).
norin.brm4 %>% summary()## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS
## Data: norin (Number of observations: 180)
## Draws: 3 chains, each with iter = 1000; warmup = 200; thin = 5;
## total post-warmup draws = 480
##
## Group-Level Effects:
## ~FISHID (Number of levels: 60)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 5.11 3.48 0.19 12.31 1.00
## sd(TRIALHypoxia) 2.46 2.22 0.08 8.34 1.00
## sd(TRIALLowSalinity) 26.24 4.50 17.65 35.41 1.00
## cor(Intercept,TRIALHypoxia) -0.11 0.50 -0.91 0.85 1.00
## cor(Intercept,TRIALLowSalinity) 0.27 0.39 -0.55 0.92 1.00
## cor(TRIALHypoxia,TRIALLowSalinity) -0.19 0.47 -0.93 0.76 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 1086 1792
## sd(TRIALHypoxia) 1831 2066
## sd(TRIALLowSalinity) 1577 1985
## cor(Intercept,TRIALHypoxia) 2120 2210
## cor(Intercept,TRIALLowSalinity) 727 812
## cor(TRIALHypoxia,TRIALLowSalinity) 623 1582
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 140.76 21.78 98.54 184.66 1.00 2262
## TRIALHypoxia -91.92 24.06 -136.82 -44.63 1.00 2368
## TRIALLowSalinity -60.10 30.53 -119.28 1.01 1.00 2079
## SMR_contr -18.13 3.65 -25.37 -10.87 1.00 2259
## MASS 0.09 0.23 -0.37 0.55 1.00 2419
## TRIALHypoxia:SMR_contr 7.99 4.60 -1.25 16.53 1.00 2287
## TRIALLowSalinity:SMR_contr 3.55 5.85 -8.03 15.03 1.00 2024
## Tail_ESS
## Intercept 2452
## TRIALHypoxia 2079
## TRIALLowSalinity 2326
## SMR_contr 2283
## MASS 2369
## TRIALHypoxia:SMR_contr 2325
## TRIALLowSalinity:SMR_contr 2290
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 18.14 1.73 14.62 21.45 1.00 1380 1553
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Conclusions:
the intercept indicates that the average metabolic change associated with the High Temperature trial (at the mean SMR_contr and MASS) across all fish was 140.76.
for the same SMR_contr and MASS, the average metabolic change was 91.92 less in the Hypoxia trial and 60.1 less in the Salinity trial than the High Temperature trial.
in the High Temperature trial, every 1 unit increase in SMR_contr was associated with a decline in metabolic change of
round(-1*norin.sum$fixed[4,1], 2)units.in the High Temperature trial (and mean SMR_contr), every 1 unit increase in MASS was associated with an increase in metabolic change of
round(-1*norin.sum$fixed[5,1], 2)units.the rate of change in metabolic change associated with the Hypoxia trial was
round(-1*norin.sum$fixed[6,1], 2)units higher than that of the High Temperature trial andround(-1*norin.sum$fixed[7,1], 2)units higher in the Low Salinity trial.the variance in intercepts across all fish is 5.11
the variance in responses between High Temperature and Hypoxia trials across all fish is 2.46.
the variance in responses between High Temperature and Low Salinity trials across all fish is 26.24.
the scale of variance between fish within a treatment
(sigma, 18.14) is substantially higher than that both within the High Temperature trial and between the High Temperature and Hypoxia trials, yet lower than that between High Temperature and Low Salinity trials.
as_draws_df (posteriors)
norin.brm4 %>% as_draws_df()